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ABSTRACT 


Hypobaric hypoxia (HH) exposure can cause serious 
brain injury as well as life-threatening cerebral 
edema in severe cases. Previous studies on the 
mechanisms of HH-induced brain injury have been 
conducted primarily using non-primate animal 
models that are genetically distant to humans, thus 
hindering the development of disease treatment. 
Here, we report that cynomolgus monkeys (Macaca 
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fascicularis) exposed to acute HH developed human- 
like HH syndrome involving severe brain injury and 
abnormal behavior. Transcriptome profiling of white 
blood cells and brain tissue from monkeys exposed 
to increasing altitude revealed the central role of the 
HIF-1 and other novel signaling pathways, such as 
the vitamin D receptor (VDR) signaling pathway, in 
co-regulating HH-induced inflammation processes. 
We also observed profound transcriptomic 
alterations in brains after exposure to acute HH, 
including the activation of angiogenesis and 
impairment of aerobic respiration and protein folding 
processes, which likely underlie the pathological 
effects of HH-induced brain injury. Administration of 
progesterone (PROG) and steroid neuroprotectant 
5a-androst-3ß,5,6ß-triol (TRIOL) significantly 
attenuated brain injuries and rescued the 
transcriptomic changes induced by acute HH. 
Functional investigation of the affected genes 
suggested that these two neuroprotectants protect 
the brain by targeting different pathways, with PROG 
enhancing erythropoiesis and TRIOL suppressing 
glutamate-induced excitotoxicity. Thus, this study 
advances our understanding of the pathology 
induced by acute HH and provides potential 
compounds for the development of neuroprotectant 
drugs for therapeutic treatment. 


Keywords: Acute hypobaric hypoxia; Cynomolgus 
monkeys; Brain injury; Neuroprotectant; Gene 
regulatory networks 


INTRODUCTION 


High altitude, which represents one of the most extreme 
environments on Earth, creates hypobaric hypoxia (HH) 
conditions to which only a small proportion of people can 
adapt. As the brain is the most hypoxia-intolerant organ, most 
people who ascend too rapidly to high altitudes cannot 
acclimatize to the HH environment and frequently suffer from 
a range of symptoms caused by acute hypoxia (Wilson et al., 
2009). In severe cases, patients can develop serious high- 
altitude cerebral edema (HACE). With over 35 million people 
around the world traveling to high-altitude regions annually for 
tourism, sport, or work (Martin & Windsor, 2008), brain 
damage caused by HH has become a health hazard for 
lowland people who move to high altitudes. So far, the main 
mechanisms underlying brain damage include metabolic 
disturbance of neural cells, increased permeability of brain 
microvasculature, and oxidative stress (Basnyat & Murdoch, 
2003; Wilson et al., 2009). However, details on the cascades 
and gene regulatory networks through which brain damage is 
induced by HH have not yet been determined. 

Several drugs, such as acetazolamide and dexamethasone, 
are proposed to treat or even prevent HH-induced brain 
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damage and edema based on increased oxygen-carrying 
capacity or anti-inflammatory effects (Imray et al., 2010). 
However, their efficacy is still controversial and numerous side 
effects have been reported in relation to their clinical use 
(Imray et al., 2010; Rabinstein, 2006). Progesterone (PROG), 
a well-known endogenous neuroprotective sex hormone, is 
reported to be associated with acute HH acclimatization 
(Stein, 2008). In addition, serum levels of PROG are found to 
increase in men who travel to high-altitude areas (Basu et al., 
1997), suggesting the possibility that PROG may ameliorate 
HH-induced brain damage. Furthermore, as first reported by 
our group (Chen et al., 2013a), 5a-androst-38,5,6B-triol 
(TRIOL) is a novel neuroprotectant and analogue of the 
endogenous neuroprotective steroid cholestane-38,5a,6B-triol 
(Hu et al., 2014), thus its effects on HH-induced brain damage 
are also worth investigation. 

Previous pathological and molecular studies of HH-induced 
brain injuries have focused primarily on small mammals such 
as mice and rats (Imray et al., 2010). However, these rodents 
display distinct hypoxia tolerance and physiological responses 
to HH in comparison to humans. For example, both acute HH 
conditions and exhaustive exercise are necessary to induce 
HACE-like symptoms in rats (Guo et al., 2013), and very rapid 
ascent speed (velocity of 50 m/s within 5 min to 6 000 m) is 
required to induce HACE-like symptoms in mice (Huang et al., 
2015). In addition, the potential drug targets in rodents may 
differ from those in humans due to the large genetic 
differences between rodents and primates after ~96 million 
years of evolutionary divergence (Nei et al., 2001). In 
particular, the central nervous system (CNS) of primates is 
much more complex than that of rodents (Lissa et al., 2013). 
Thus, given the limitations of rodent models, the establishment 
of non-human primate models for HH study is urgently needed 
and should provide valuable insights into the molecular 
cascades and gene regulatory networks underlying HH- 
induced brain diseases and further our understanding of the 
effects and molecular mechanisms of neuroprotectants as 
potential drugs. 

In the current study, we elucidated the spatial and temporal 
influence of acute HH on gene expression and examined the 
possible effects of PROG and TRIOL on HH-induced brain 
damage in vivo. Specifically, we established a non-human 
primate model of HH-induced brain damage and profiled the 
transcriptomes of white blood cells (WBCs) and brain tissue 
from cynomolgus monkeys before and after acute HH 
exposure (i.e., increasing altitude). Extensive gene regulatory 
analyses revealed a dynamic change in the WBC 
transcriptome response to HH as well as the gene regulatory 
network of newly identified transcriptomic hub genes. 
Moreover, the application of the two neuroprotectants 
effectively protected the brains of cynomolgus monkeys from 
HH-induced injury. Finally, PROG and TRIOL exerted their 
effects via different pathways, the former through 
erythropoiesis and the latter by suppression of glutamate- 
induced excitotoxity. 


MATERIALS AND METHODS 


Experimental animals 

All experiments were conducted in accordance with the 
Chinese Laws for the Protection of Animals. The experimental 
protocols were approved by the Ethics Committee of 
Zhongshan School of Medicine, Sun Yat-Sen University 
according to the ARRIVE (Animal Research: Reporting of /n 
Vivo Experiments) guidelines (Kilkenny et al., 2010). All 
animal-based procedures were performed in strict adherence 
to the National Standards of Treating Experimental Animals 
(2006 version). Efforts were taken to minimize suffering and to 
ensure the welfare of monkeys during experimentation. 

Twenty-four male cynomolgus monkeys (6.0-—6.5 years old, 
6.8-7.5 kg) were obtained from Gaoyao Kangda Laboratory 
Animals Science & Technology Co., Ltd, Zhaoqing, 
Guangdong Province, China. The animals were transported to 
the Third Military Medical University, Chongqing, China. 
Monkeys were singly housed in cages in a controlled 
environment at a temperature of 25+1 °C, relative humidity of 
60%, and circadian 12 h light/dark cycle and were fed 
routinely. 

For experimentation, monkeys were randomly divided into 
four groups, each containing six individuals. The first group 
was maintained under normobaric normoxia (NN) conditions 
(at an altitude of 320 m) as a control. The other three groups 
were subjected to hypobaric hypoxia (HH), as detailed below. 
Groups 3 and 4 were treated with progesterone (HH+PROG) 
and 5a-androst-38,5a,6B-triol (HH+TRIOL), respectively. 

After 48h at a simulated altitude of 7500m, the HH 
monkeys were sacrificed. The NN monkeys were also 
sacrificed after 3 d in the NN environment. All monkeys were 
euthanized by bloodletting from the carotid under anesthesia 
using a mixture of injectable ketamine hydrochloride 
(0.06 mg/kg) and xylazine hydrochloride (0.02 mg/kg). 


Pharmacological treatments 

Intravenous (10 mg/mL) and extended-release intramuscular 
injections (50 mg/mL) of TRIOL were provided by Guangzhou 
Cellprotek Pharmaceutical Company Co., Ltd. (China). The 
TRIOL (Batch No. 101124) was dissolved in vehicle containing 
0.9% sodium chloride and 20% hydroxypropyl-B-cyclodextrin 
(HP-B-CD) for intravenous injection, and in vehicle containing 
12% glycerin, 20% HP-B-CD, and 0.19% CMC-Na_ for 
extended-release intramuscular injection. Injectable 
progesterone (20 mg/mL) was purchased from Zhejiang 
Xianju Pharmaceutical Co., Ltd. (China), with a CFDA 
ratification No. of Guo Yao Zhun Zi-H33020828. The HP-B-CD 
(Batch No. 100309) was purchased from Xi'an Deli Biology & 
Chemical Industry Co., Ltd. (China). Injectable ketamine 
hydrochloride (150 mg/mL) was purchased from Shenyang 
Veterinary Drugs Co., Ltd. (China) with an approval No. of 
(2011)060022668 by the Ministry of Agriculture of the People's 
Republic of China (PRC). Injectable xylazine hydrochloride 
(100 mg/mL) was purchased from Dunhua_ Shengda 
Veterinary Drugs Co., Ltd. (China) with an approval No. of 
(2010)070031581 by the Ministry of Agriculture of the PRC. 


Acute hypobaric hypoxia and neuroprotective steroid 
agent treatments 

For HH treatment, the simulated altitude in the hypobaric 
chamber was increased from 320 m to 6 000 m at a velocity of 
3 m/s, and from 6 000 m to 7 500m at a velocity of 2 m/s. 
Altitude simulation was started at 320m (stage A), then 
suspended at 3000 m (stage B), 4 500m (stage C), and 
6 000 m (stage D) for 50 min each, and finally maintained at 
7 500 m for an initial 24 h (stage E) and then another 24h 
(stage F) for a total of 48h. At every stage, monkeys were 
maintained under HH conditions for 30min at the 
corresponding altitudes, and were then subjected to blood 
sample collection and drug administration within 20 min. 
Because of safety concerns for the experimenters, the altitude 
was decreased from 7 500 m to 6 000 m for blood collection 
and drug administration at stage E within 30 min. Throughout 
the HH experimental exposure, the temperature, relative 
humidity, and air-flow velocity of the chamber were maintained 
at 22 °C, 60%, and 5 L/min, respectively. 

For the HH+PROG group, animals were intramuscularly 
injected with 15 mg/kg of PROG three times, first at 12 h 
before HH treatment, then during stages D and E. For the 
HH+ TRIOL group, animals were intravenously injected with 
10 mg/kg of TRIOL 5 min before HH treatment, and at stages 
B, C, D, and E. Owing to TRIOL's short half-life (t,;2) of 0.5 h, 
the extended-release intramuscular injection of TRIOL 
(30 mg/kg) was also administered at stages D and E. 


Assessment of skeletal muscle coordination behavior 
Behavioral assessment was performed by trained 
experimenters blind to group information. Monkey behavior 
was monitored 2h after arrival at the simulated altitude of 
7 500 m. We observed their activities over the next 30 min and 
employed a modified scale method to assess skeletal muscle 
coordination, with higher scores indicating more severely 
impaired behavioral states. The scoring system for monkey 
behavior assessment is presented in Supplementary Table 
S1. 


Measurement of brain water content 

Left brain hemispheres were used to measure brain water 
content with a precision electronic scale (Sartorius, BSA224S, 
Germany). Briefly, whole left hemispheres were immediately 
removed and measured as 'wet weight'. The left hemispheres 
were then placed in an electric thermostatic oven at 60 °C and 
weighed repeatedly every day until reaching a constant value 
(dry weight). The percentage of brain water content was 
calculated as: water content (%)=(wet weight-dry weight)/wet 
weightx100%. 


Transmission electron microscopy (TEM) 

The frontal cortices of brains were divided into 1 mmx1 mmx 
1 mm pieces and fixed in 2.5% glutaraldehyde overnight at 
4°C. The fixed brain tissues were dehydrated through an 
ethanol series, embedded in epoxy resin, and cut into 60 nm 
ultrathin sections with a Leica EM UC6 Ultramicrocut 
(Germany). The sections were mounted on copper grids and 
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then stained in uranyl acetate and citric acid. Images were 
captured using a transmission electron microscope (Tecnai 
G2 Spirit Twin, FEI, USA). The experimental procedures for 
TEM were completed in a double-blind manner. 


Nissl staining 

The frontal cortices of brains were cut into 5 mm-thick pieces 
immediately after removal and then fixed in 4% 
paraformaldehyde. All samples were dehydrated in an 
increasing series of ethanol, cleared in xylene, and embedded 
in paraffin. Brain sections of 5 um-thickness were cut using a 
rotary microtome (Shandon Finesse 325, Thermofisher 
Scientific, USA) and then deparaffinized with xylene and 
hydrated, followed by Nissl staining (0.5% cresyl violet) for 
histopathological assessment of neuronal injury. For Nissl 
staining, injured neurons were characterized by dark staining, 
condensed nuclei, shrunken cell bodies, or weak staining with 
irregular shapes. Images were acquired with a bright-field 
microscope (Nikon ECLIPSE Ti-U, Japan). 


WBC isolation 

WBCs were isolated from whole blood samples using Red Cell 
Lysis Buffer (Tiangen, RT122, China). In brief, red blood cells 
(RBCs) were lysed by gently mixing the blood sample with 
Red Cell Lysis Buffer (1:3), followed by incubation on ice for 
5 min. The precipitate was then harvested by centrifugation at 
2 500 g for 5 min at 4 °C. The WBC samples were ready for 
total RNA extraction after the isolation procedure was 
repeated three times. 


RNA extraction, library construction, and sequencing 
Total RNA was extracted and purified from WBCs and frontal 
cortices using Trizol reagent (Invitrogen, 15596, USA) 
according to the manufacturer's standard protocols. RNA 
quality and concentration were assessed using an Agilent 
Bioanalyzer 2100 and RNA 6000 Nano Kit (Agilent 
Technologies, USA). One of the TRIOL-treated monkeys was 
removed from RNA-seq analysis because the blood and brain 
RNA samples were poorly preserved due to depletion of dry 
ice during delivery to the sequencing company. One of the 
stage-A WBC samples from the HH+PROG group was also 
removed due to severe degradation of RNA. 

For library construction, oligo-dT-coupled beads were first 
used to enrich poly-A+RNA molecules. First-strand cDNA 
synthesis was performed using random hexamers and 
Superscript Il reverse transcriptase (Invitrogen, USA). 
Second-strand cDNA synthesis was performed using E. coli 
DNA Poll (Invitrogen, USA). A Qiaquick PCR purification kit 
(Qiagen, Germantown, MD, Germany) was used to purify the 
double-strand cDNA. cDNA was sheared with a nebulizer 
(Invitrogen, USA) into 100-500 bp fragments. The fragments 
were ligated to a Lumina PE adapter oligo mix after end repair 
and the addition of a 3' dA overhang. The 310-350 bp 
fragments were then collected by gel purification. After 15 
cycles of PCR amplification, the libraries were subjected to 
paired-end sequencing using an HiSeq2500 (Illumina, USA) 
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platform. Finally, 357.58 Gb of transcriptome sequences were 
generated from the 91 samples for transcriptome analysis. 


Measurement of intraneuronal calcium concentration 
([Ca**Ji) 

The intracellular [Ca?*]i of cultured primary rat cortical neurons 
was measured as per Minta et al. (1989). Briefly, after 8 d of in 
vitro culture on glass cover slips, primary cortical neurons 
were washed once with Hank's balanced saline solution 
(HBSS) loading buffer. We added 100 uL of HBSS loading 
buffer containing 5 pmol/L Fluo-3/AM (Molecular Probes, 
F1241, USA) dropwise onto the cover slips, which were then 
incubated with neurons for 45 min at room temperature before 
being washed twice with HBSS testing buffer to remove 
residual fluorescent dye. Neurons were then transferred to the 
perfusion chamber of an OctaFlow Perfusion System (ALA 
Scientific Instruments, Farmingdale, NY, USA), which was 
placed on the stage of a TCS SP2 laser scanning confocal 
microscope (Leica Microsystems, Mannheim, Germany). 
Using an excitation wavelength of 488 nm and emission 
wavelength of 526 nm, at least 10 randomly selected neurons 
on every slide were scanned continuously at 6 s intervals to 
determine real-time fluorescence intensity. Dynamic changes 
in the fluorescence intensity of every neuron were recorded 
automatically by the Leica confocal system. Changes in [Ca?*Ji 
(A[Ca?*]i) were represented by changes in fluorescence 
intensity before and after drug administration and were 
expressed as: A[Ca**]i=(F (peak fluorescence intensity after 
administration) — FO (baseline fluorescence intensity before 
administration))/FOx100%. 


Statistical analysis of behavioral 
experiments 

Statistical analyses were performed using SPSS 19.0 and 
graphs were drawn using GraphPad Prism 6.0. Most data, 
including brain water content and number of injured neurons, 
were analyzed using Student's t-test or one-way analysis of 
variance (ANOVA) with Dunnett's post-hoc test for multiple 
comparisons. Behavioral impairment scores were analyzed 
using the Kruskal-Wallis test. All data are represented as 
means+SD unless otherwise noted. P-values of less than 0.05 
were considered statistically significant. 


and neurological 


Source of genome reference and gene annotation 
The reference genome and gene annotation of the 
cynomolgus monkey used in this study were originally 
generated by Yan et al. (2011) and downloaded from GigaDB 
(http://gigadb.org/dataset/100003). Specifically, the following 
files were downloaded: 

Genome assembly: ftp://parrot.genomics.cn/gigadb/pub/ 
10.5524/100001_101000/100003/CE.cns.all.fa.gz 

Gene annotation (GFF): ftp://parrot.genomics.cn/gigadb/pub/ 
10.5524/100001_101000/100003/CE. gff.gz 

Gene annotation (CDS): ftp://parrot.genomics.cn/gigadb/ 
pub/10.5524/100001_ 101000/100003/CE.cds.fa.gz 

Gene annotation (PEP): ftp://parrot.genomics.cn/gigadb/pub/ 
10.5524/100001_101000/100003/CE. pep.fa.gz 


Transcriptome 
measurement 
The Hisat2 (v2.0.4) (Kim et al., 2015) package was employed 
to align transcriptome reads to the reference genome of 
cynomolgus monkey with default parameters. Gene 
expression was measured as RPKM (reads per kilobase of 
transcript per million mapped reads) using an in-house script 
(supplementary file "count_reads_num. pl"), and library sizes 
were normalized with the geometric method in DESeq2 
(v1.10.1) (Love et al., 2014) (Supplementary Data). Genes 
expressed robustly (i.e., mean RPKM>1) in at least one 
sample group were used for principal component analysis 
(PCA). 


alignment and expression level 


Identification of differentially expressed genes (DEGs) 

We applied a strict protocol to identify DEGs between groups 
of samples. Candidate DEGs were first identified using 
DESeq2 (v1.10.1) (Love et al., 2014), edgeR (v3.12.1) 
(Robinson et al., 2010), and Cuffdiff (v2.2.1) (Trapnell et al., 
2013) separately. For DESeq2, the default method was used 
to normalize library sizes, the fitting type of dispersions to 
mean intensity was set as "parametric", and the method to test 
DEG significance was set as "nbinomWaldTest". For edgeR, 
the trimmed mean of M-values (TMM) normalization method 
was used to normalize library sizes, the quantile-adjusted 
conditional maximum likelihood (CML) method was used to 
estimate dispersions, and the genewise exact test 
("exactTest") was used to calculate the significance of DEGs. 
For Cuffdiff, the method for library size normalization was set 
as "--geometric-norm" and other parameters were set as 
default. P-values were adjusted to allow for multiple testing by 
the Benjamini-Hochberg False Discovery Rate (FDR) 
(Benjamini & Hochberg, 1995). A DEG was required to meet 
the following criteria: (1) at least two of the three DEG- 
detection software reported a FDR<0.05; (2) the lower quartile 
of expression level in the up-regulated group was greater than 
the upper quartile of expression level in the down-regulated 
group; and, (3) the mean RPKM was not less than 5 in the up- 
regulated group. Of note, for convenience, the FDR values 
shown for the genes mentioned in the Results section were all 
from edgeR. 

To verify the robustness of our DEG-detection strategy and 
to calculate the false-positive rate, samples were randomly 
divided into pseudo groups. For WBC samples, three pseudo 
groups were generated, each consisting of two samples 
randomly chosen from stage A, stage D, and stage F, 
respectively. For brain samples, two pseudo groups were 
generated, each consisting of three samples randomly chosen 
from the HH and NN groups, respectively. The procedures for 
DEG detection described above were performed, and the 
numbers of DEGs between pseudo groups of the same tissue 
were considered as false positives. In contrast to the 
thousands of significant DEGs identified between real groups, 
no DEGs were identified between pseudo groups for WBCs 
and brain tissue, confirming the robustness of our strategy and 
reliability of our results. 


Gene ontology annotation, Kyoto Encyclopedia of Genes 
and Genomes (KEGG) annotation, and enrichment 
analysis 
Gene ontology (GO) terms (Ashburner et al., 2000) for the 
cynomolgus monkeys were annotated according to their 
orthologous relationships with rhesus macaques and humans 
(Supplementary Table S5). KEGG pathways (Kanehisa & 
Goto, 2000) were annotated employing online KAAS (KEGG 
Automatic Annotation Server) tools (www.genome.jp/ 
tools/kaas/) (Moriya et al., 2007) using all protein sequences 
of the cynomolgus monkeys as inputs with a BBH (bi- 
directional best hit) method (Supplementary Table S5). 
Fisher's exact tests were employed to identify whether a list 
of genes (foreground genes) was enriched in specific GO 
terms or KEGG pathways, with comparisons of the number of 
foreground genes annotated to the specific GO/KEGG, 
number of foreground genes not annotated to the specific 
GO/KEGG, number of background genes (excluding 
foreground genes) annotated to the specific GO/KEGG, and 
number of background genes (excluding foreground genes) 
not annotated to the specific GO/KEGG. P-values were 
adjusted to allow for multiple testing by means of the 
Benjamini-Hochberg False Discovery Rate (Benjamini & 
Hochberg, 1995). 


Co-expression module and hub gene identification 

The WGCNA (v1.51) (Langfelder & Horvath, 2008) package in 
R (v3.2.3) was employed for co-expression network analyses 
using the percentage of transformed RPKM of all robustly 
expressed genes. Co-expression modules were identified 
using standard methods, with correlation between genes 
calculated by the Pearson method, and the distance for 
hierarchical clustering calculated by the Ward method (Ward 
Jr, 1963). Similar modules were merged by WGCNA with the 
parameter MEDissThres=0.25. The eigengene (vector 
presenting expression dynamics of a module between 
samples) and kME (correlation of the gene with corresponding 
module eigengene) in every module were also calculated by 
WGCNA. For every module, genes with kME>0.9 were 
considered hub genes. 


Analysis of impact of drugs on transcriptome dynamics 
For WBCs, the sequencing of the HH and drug-treated groups 
was not performed in the same batch. Hence, laboratory 
conditions, personnel differences, and reagent lots could have 
caused variations between the groups, i. e., batch effects 
(Leek et al., 2010). It should be noted that stage A samples 
from the PROG and TRIOL groups were collected prior to 
drug injection; therefore, the stage A DEGs between the HH 
and drug-treated groups were likely a reflection of batch 
effects rather than drug-specific effects, and were thus 
excluded from subsequent analyses on DEGs between HH 
and drug-treated groups in stages D and F. 

For brains, we identified the DEGs between the NN, HH, 
and drug-treated groups in pairwise comparisons. Firstly, we 
focused on the impact of drug treatment on HH-responding 
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genes, which were DEGs between the NN and HH groups. If 
an HH-responding gene was a DEG between the HH and 
drug-treated groups but not a DEG between the NN and drug- 
treated groups (i. e., the expression level of the gene in the 
drug-treated group was not significantly different from that in 
the NN group but was significantly different from that in the HH 
group), the gene was considered to be a drug strongly 
responding gene (SRG). If an HH-responding gene was a 
DEG between the NN and drug-treated groups but not a DEG 
between the HH and drug-treated groups (i.e., the expression 
level of the gene in the drug-treated group was not 
significantly different from that in the HH group but was 
significantly different from that in the NN group), the gene was 
considered to be a drug non-responding gene (NRG). The 
other HH-responding genes with expression levels that fell 
between the drug-treated, NN, and HH groups were 
considered drug weakly responding genes (WRGs). If a gene 
was not significantly differentially expressed by HH in contrast 
to NN (i.e., not a DEG between the HH and NN groups) but 
was significantly differentially regulated by a drug in contrast 
to both NN and HH conditions, it was considered to be a drug- 
induced DEG. Even if a gene was significantly up- or down- 
regulated by HH, if that gene was further significantly up- or 
down-regulated but not restored by the drug, it was also 
considered to be a drug-induced DEG, not an NRG. 


RESULTS 


Monkeys develop severe brain injuries under acute HH 
conditions that can be significantly alleviated by PROG 
and TRIOL 

Twenty-four healthy, adult male cynomolgus monkeys were 
classified into four groups (n=6): i. e., normobaric normoxia 
(NN), HH, HH+PROG, and HH+TRIOL. To characterize the 
behavioral and pathological changes caused by acute HH and 
to assess the drug effects, animals from the HH, HH+PROG, 
and HH+TRIOL groups were placed in a hypobaric chamber, 
which mimicked ascending altitudes from 320 to 7 500 m ata 
velocity of 2-3 m/s (HH group). Altitude simulation for each 
monkey started at 320 m (stage A) and was then sequentially 
suspended at 3 000 m (stage B), 4 500 m (stage C), and 
6000m (stage D) for 50min each. Altitude was then 
maintained at 7 500 m for 24 h (stage E) and then a further 24 h 
(stage F) for a total of 48 h (Figure 1A). During HH treatment, 
monkeys in the HH+PROG and HH+TRIOL groups were given 
either PROG or TRIOL, respectively (Figure 1B, C; see 
Methods). In parallel, monkeys in the NN group were placed in 
a normobaric normoxia environment at an altitude of 320 m to 
serve as controls. 

The impact of acute HH exposure on behavior was 
assessed 2 h after ascending to 7 500 m. In sharp contrast to 
the NN monkeys, which exhibited high activity levels with 
frequent walking, climbing, and eating, the HH monkeys lost 
their balance and lay prostrate with limited body movements, 
or even remained in lateral or dorsal recumbency during the 
entire behavioral assessment. All six monkeys in the HH 
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group displayed anorexia, vomiting, motor deficits, and ataxia. 
In contrast, monkeys in the HH+PROG and HH+TRIOL 
groups were able to maintain their balance and remain seated, 
with frequent head and forelimb movements. Several PROG 
or TRIOL-treated animals exhibited normal feeding behavior 
or short-distance movements (a few steps), although vomiting 
still occurred in some monkeys (two PROG-treated monkeys 
and four TRIOL-treated monkeys). We then employed 
behavioral deficit scoring to obtain a quantitative measure of 
skeletal muscle coordination by recording the movements of 
each monkey over 30 min (Supplementary Table S1) (Kito et 
al., 2001). All NN monkeys scored 0 (i. e., walked normally), 
whereas the HH monkeys scored between 14 (i. e., could not 
stand or sit for long) and 18 (i. e., no movement). PROG and 
TRIOL-treated HH monkeys scored between 10 (i. e., able to 
stand spontaneously) and 14; this was significantly lower than 
the scores of untreated HH monkeys (ANOVA P<0.05, Figure 
1D), indicating that PROG and TRIOL were able to partially 
restore motor coordination in HH monkeys. 

Cerebral edema is normally manifested as the increase of 
the water content in the brain tissue. We found that brain 
water content was significantly higher in the HH group 
(76.71%+0.30%) than in the NN group (76.16%+0.26%) 
(ANOVA P<0.01, Figure 1E). The proportion of brain water 
content was significantly reduced in both the HH+PROG and 
HH+TRIOL groups to 76.12%+0.46% and 76.28%+0.23%, 
respectively, similar to levels in the NN group (Figure 1E). 
Additionally, the pericapillary space was markedly wider in 
untreated HH brains than in other groups (Figure 1F), 
suggesting that blood-brain barrier (BBB) disruption and 
vasogenic edema, which are considered to cause cerebral 
edema in patients with HACE (Bartsch & Swenson, 2013), had 
been induced by HH and reversed by drug administration. 
Nissl staining of the frontal cortex revealed many injured 
neurons in the HH group only. These cells were characterized 
by dark staining, condensed nuclei, shrunken cell bodies, or 
weak staining with irregular shapes (Figure 1G). The 
proportion of injured neurons in the HH group was 21.1%, 
significantly higher than that in the NN group (1.4%) (ANOVA 
P<10-). This ratio was restored to 7.5% in the HH+PROG 
group (ANOVA P<0.01) and 6.7% in the HH+TRIOL group 
(ANOVA P<0.01, Figure 1H; see Methods). 

Taken together, these findings suggest that cynomolgus 
monkeys develop severe motor impairments and brain 
damage after acute exposure to HH environments, and this 
can be significantly mitigated by either PROG or TRIOL. 


WBC dynamic transcriptome analysis reveals three major 
regulatory modules responding to acute HH 

To systematically investigate the genetic regulatory networks 
that participate in the HH response, we sequenced and 
analyzed the transcriptome of WBCs obtained at each of the 
six HH stages, A-F (n=6, 36 WBC samples in total). Based on 
PCA of 12 198 robustly expressed genes, the degree of HH 
exposure (i. e., a function of altitude plus exposure time), 
rather than individual differences, appeared to be the main 
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Figure 1 Acute HH-induced behavioral and cerebral impairments were attenuated by PROG and TRIOL treatments 
A: Experimental procedure for HH treatment of cynomolgus monkeys (HH group, n=6). Vertical bars along y-axis indicate medical definitions of 
high, very high, and extreme altitudes. Horizontal bars represent duration spent at each altitude. Blood drop and brain icons indicate altitudes and 
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time points at which blood and brain samples, respectively, were collected for RNA-seq analyses. B: Experimental procedure for drug treatment 
(HH+ PROG and HH+TRIOL groups, each n=6). Needle tubes indicate altitudes and time points at which drugs were injected; blood drop and brain 
icons indicate altitudes and time points at which blood and brain samples were collected, respectively. Drugs were injected after blood collection at 
each stage. Samples from stages highlighted in red were selected for RNA-seq analyses. C: Chemical structures of two neuroprotectants (PROG 
and TRIOL) used in this study. D: Mean skeletal muscle coordination scores for each experimental group, with higher scores representing a lower 
degree of coordination, assessed after 2 h at 7 500 m. Each group was compared with HH group (*: Kruskal-Wallis test P<0.05). Data are presented 
with error bars indicating means+SD, n=6 per group. E: Mean brain water content for each experimental group, measured as 1—dry/wet weight of 
left hemisphere. Each group was compared with HH group (*: one-way ANOVA P<0.05; ™: P<0.01, HH n=5, each other group n=6). Error bars 
represent SD. F: Representative images showing ultrastructure of capillaries in frontal cortex of each experimental group, taken after 48 h at 7 500 
m. Blood-brain barrier (BBB) disruption and vasogenic edema evident in HH group but not in other groups, as characterized by severely shrunken 
capillaries with fluid penetrating into pericapillary space. Red arrowheads indicate capillary boundary. Scale bar: 2 um. G: Representative images 
showing Nissl staining of frontal cortex tissues of each experimental group. Neuronal injuries evident in HH group and markedly attenuated in drug- 
treated groups. Black squares in first row denote areas magnified in second row. Red arrowheads indicate injured cells with dark staining, 
condensed nuclei, shrunken cell bodies, or weak staining with irregular shapes. Red scale bar: 75 um, black scale bar: 20 um. H: Percentage of 
injured neurons in frontal cortex of each experimental group based on Nissl staining shown in panel G. Each group was compared with HH group (”: 


one-way ANOVA Ps0.01; “”: Ps0.001, NN n=5, all other groups n=4). Error bars represent SD. 


factor affecting the genome-wide WBC gene expression 
patterns (Figure 2A). We then identified DEGs based on 
comparisons between any two of the six HH stages (see 
Methods) and obtained a total of 5 174 DEGs across all 
comparisons (ranging from 0-3 864 for different comparisons). 
As stage A (simulated altitude of 320 m) represents the NN 
condition, the number of DEGs between this and other stages 
reflects the degree of gene expression change in response to 
varying levels of HH exposure. The number of DEGs 
increased incrementally from stages B to D, reaching its 
highest at stage D. During stages E and F, the number 
decreased (Figure 2B), implying that the increased gene 
expression changes associated with HH were reversed in later 
stages, probably because of organismal acclimatization to 
hypoxia after 24 h of exposure to HH (Berglund, 1992). 
Co-expression analysis was next used to identify co- 
expression modules associated with HH changes (see 
Methods). Eight modules were identified, ranging in size from 
64 to 5 607 genes (Supplementary Figure S1A, B). Three 
modules were significantly enriched with DEGs (Fisher's exact 
test P<10-'5) and showed strong correlation between the 
expression pattern and progression of HH (ANOVA P<107; 
Supplementary Table S2). Hence, these were considered the 
primary HH-responding modules (HH modules). The first HH 
module contained 3 114 genes, 2 123 (68%) of which were 
DEGs (Supplementary Table S2). This module (termed the 
‘impulse module’) showed an ‘impulse pattern’, whereby the 
expression levels were monotonically increased or decreased 
in the early stages of HH (stages A-D), but recovered to a 
certain extent during the later stages (stages E-F) (Figure 
2C). The second HH-responding module (termed the 
‘sustained module’) contained 571 genes, 479 (84%) of which 
were DEGs (Supplementary Table S2). In this module, 
changes in expression levels of genes that responded to 
early-stage changes in altitude were sustained in the later 
stages of HH (Figure 2C). The third module (termed the ‘late- 
responding module’) contained 552 genes, 398 (72%) of 
which were DEGs (Supplementary Table S2). Genes in this 
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module showed no differential expression until the later stages 
of HH (Figure 2C). 


Impulse responses of hypoxia/HIF-1 mediated innate 
immunity and inflammatory processes to acute HH 

Genes in the impulse module were significantly implicated in 
several classical pathways critical to hypoxia response and 
inflammation (Figure 2D-E, and Supplementary Figure S1C). 
As noted in previous small mammal and human cell line 
studies, we confirmed that the HIF-1 signaling pathway also 
played a regulatory role in adaptive responses promoting 
optimal oxygen utilization (Semenza, 2009) and in 
inflammatory responses (Palazon et al., 2014) to hypoxia in in 
vivo primate experiments. Other innate immune and 
inflammatory pathways such as the NF-kB, Toll-like receptor, 
and FoxO signaling pathways, were also activated, probably 
forming a positive feedback loop with HIF-1 (Han et al., 2016). 
The inflammatory factor IL1B and chemokine CXCL1, which 
constitute links with HIF-1 and NF-kB signaling (Hatfield et al., 
2010; Jung et al., 2003), were up-regulated by HH in the 
impulse module and are known to contribute to vascular 
impairment and vasogenic edema (Stamatovic et al., 2006), 
as observed in the brains of the HH-treated group (Figure 1F). 
Collectively, HIF-1/hypoxia pathways mediating the 
inflammatory response may contribute to pathogenesis in 
monkeys and were expressed in the impulse module during 
HH progression (Figure 2E). In contrast, the sustained module 
included genes mostly involved in adaptive immune responses 
(Figure 2D and Supplementary Figure S1C). Specifically, most 
key T cell markers, T cell receptors, and their specific kinases 
were consistently down-regulated by HH, suggesting overall 
suppression of the adaptive immune system in the HH 
environment. 


Emergence of RBC-associated inflammatory processes in 
late stages 

The late-responding module mainly contained genes involved 
in erythropoiesis and platelet activation (Figure 2D and 
Supplementary Figure S1C). Several genes encoding key 
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Figure 2 Transcriptomic dynamics of WBCs in response to acute HH. 

A: Principal component analysis (PCA) using genes robustly expressed in at least one of six HH stages (A-F). Colors denote stages of HH at which 
WBC samples were collected, as described in Figure 1A. B: Numbers of differentially expressed genes (DEGs) between stage A and each of the 
other five stages. C: Expression patterns of DEGs in three HH-responding modules, shown as log2-transformed reads per kilobase of transcript per 
million mapped reads (RPKM) fold changes of DEGs in each stage compared to stage A. Red and light blue lines indicate expression patterns of 
individual DEGs up- and down-regulated by HH, respectively. Bold lines show mean expression levels of all up- or down-regulated DEGs, with error 
bars indicating SD. UP and DOWN indicate exact number of up- or down-regulated DEGs by HH in each module. D: Kyoto Encyclopedia of Genes 
and Genomes (KEGG) pathways enriched in three HH-responding modules. Only significantly enriched pathways with a false discovery rate (FDR)- 
adjusted P-value of <0.05 were plotted. Dot color denotes P-value after FDR correction, and dot size denotes ratio of DEGs versus all expressed 
genes in each pathway. Different background colors denote classification of pathways. E: Interaction of pathways and symptoms according to 
reported articles. Blue color denotes pathways of impulse module; red color denotes pathways of late-responding module. Solid arrow line denotes 
reported direct interaction; dashed arrow line denotes reported associated relationship. F: Local co-expression network revealing genes in vitamin D 
receptor (VDR) complex as hub genes. Dark blue circles denote DEGs in impulse module; red circles denote genes in VDR complex and HIF1A as 
hub genes. 
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factors in erythropoiesis, including erythroid transcription 
factors (e.g., GATA1 and KLF1), erythropoietin receptor 
EPOR, alpha and beta-globin genes (e.g., HBA1, HBA2, and 
HBB), hemoglobin stabilizing protein AHSP, key synthases of 
heme (e.g., HMBS, ALAD, and UROS), tropomodulin-1 
TMOD1, and E3 ubiquitin-protein ligase TRIM58, were 
significantly up-regulated during stages E and F 
(Supplementary Figure S1D). This suggests that the 
enhanced expression of erythropoiesis genes was induced as 
an adaptive response after 28—52 h (stages E-F, Figure 1A) 
of HH exposure, which can regulate erythropoiesis and 
improve oxygen supply, thereby attenuating the impact of 
hypoxia in later stages. Key genes of blood coagulation and 
inflammation were also significantly up-regulated in late 
stages, including GP1BA, F2R, and other important 
procoagulant/antifibrinolytic factors (e.g., VWF, SERPINE7, 
and F13A1) that are associated with impairments of RBCs and 
vascular inflammation (Bergmeier et al., 2006; Gemmati et al., 
2016; Sparkenbaugh & Pawlinski, 2013). In addition, the 
expression of the haptoglobin gene HP, a major scavenger of 
free hemoglobin, was also significantly up-regulated by HH 
(Supplementary Figure S1D), reflecting a homeostatic 
response to the increased release of free hemoglobin from 
impaired RBCs. 

Intriguingly, in the late-responding module, we also found 
genes encoding critical markers of neutrophil degranulation, 
including neutrophil elastase (ELANE), myeloperoxidase 
(MPO), and cathepsin G (CTSG), the expression levels of 
which were up-regulated by HH (Supplementary Figure S1D). 
These activated neutrophils are tightly connected with BBB 
disruption, vascular inflammation, and vasogenic edema in 
injured brains (Ikegame et al., 2010; Liu et al., 2018), which 
could contribute to the development of the HACE-like 
symptoms observed in monkeys. 


Vitamin D receptor (VDR) signaling is a novel key 
regulator involved in HH response 

We next searched for intra-modular hub genes (i. e., genes 
with the highest degree of connectivity with other genes in the 
same module) in the three HH-responding modules to identify 
core regulators of the regulatory networks (see Methods). A 
total of 326 hub genes were identified (Supplementary Table 
S3), including those encoding key transcription factors, 
kinases, and receptors from the HIF-1 and FoxO signaling 
pathways (e.g., STAT3, FOXO4, and PTEN in the impulse 
module), T cell receptor signaling (e.g., LCK and PTPN4 in the 
sustained module), and erythrocyte development and 
maintenance (e.g., TMOD1, TRIM58, and KLF7 in the late- 
responding module). 

Interestingly, our hub gene analysis also identified the 
vitamin D receptor (VDR) signaling pathway as an important 
but previously unreported pathway involved in the HH 
response. Specifically, genes encoding two key transcription 
factors, VDR and retinoid X receptor alpha (RXRA), and two 
co-activators, nuclear receptor coactivator 1 (NCOA1) and 
tripartite motif containing 24 (TRIM24), were identified as hub 
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genes in the impulse module and were up-regulated at an 
early stage (Figure 2F). VDR and RXRA encode proteins that 
together form a heterodimeric complex that serves as a 
nuclear receptor for calcitriol, the hormonally active metabolite 
of vitamin D (Kato, 2000). Previous studies have revealed that 
VDR can inhibit the transcriptional activity of NF-KB (Chen et 
al., 2013b) and HIF-1 (Chung et al., 2009). VDR may also be 
up-regulated by Toll-like receptors and participate in the 
regulation of both the innate and adaptive immune systems 
(Adams & Hewison, 2008). These observations raise the 
possibility that the VDR pathway may function as a hub 
regulatory network that balances HIF-1 signaling and its 
inflammatory responses. Consistent with this, the vitamin D- 
binding protein in humans is associated with HH adaption in 
high-altitude dwelling native peoples (Ahmad et al., 2013), 
whereas VDR in mouse endothelial cells protects brains from 
hypoxia/reoxygenation-induced BBB disruption (Won et al., 
2015). 


Diverse effects of PROG and TRIOL on HH-induced WBC 
transcriptome dynamics 

To evaluate the effects of two neuroprotective steroid agents 
on the blood cell transcriptome, we collected WBC samples 
for RNA-seq analysis from the HH+PROG and HH+TRIOL 
groups at stages D and F (Figure 1B), because stage D 
exhibited the most dynamic response to HH (Figure 2B) and 
stage F was the final stage of our analysis. WBC samples 
were also collected from the same animals at stage A, before 
drug injection, to serve as controls (Figure 1B). PCA revealed 
that neither PROG nor TRIOL treatment of HH animals had a 
significant impact on their WBC transcriptome dynamics, and 
the degree of HH exposure remained the primary 
differentiating factor between samples (Figure 3A). 

To further quantify the potential effects of the two 
neuroprotective agents on specific genes and pathways, we 
identified agent-induced DEGs by comparing the gene 
expression between samples from agent-treated (HH+PROG 
or HH+TRIOL) and untreated HH groups at stages D and F. 
We identified 134 DEGs induced by PROG in stages D (70 
DEGs) and F (72 DEGs) and found that these DEGs 
overlapped significantly with those from the late-responding 
module (Figure 3B and Supplementary Figure S2A, Fisher's 
exact test: P<101?). In particular, key erythrocyte-associated 
genes were among the DEGs with the largest expression 
changes (top DEGs) in stage F and significantly up-regulated 
by PROG, including the alpha and beta-globin genes (HBAT7: 
FC=3.03, FDR=0.000 5; HBA2: FC=2.99, FDR=0.000 8; HBB: 
FC=3.40, FDR=0.0001), heme synthases (SLC25A39: 
FC=2.10, FDR=0.008 3; ALAS2: FC=3.43, FDR=0.000 3), and 
bisphosphoglycerate mutase (BPGM: FC=2.33, FDR=0.026 6) 
(Figure 3C). Thus, although PROG appeared to have a 
relatively subtle effect on the WBC transcriptomes, we 
speculate that this steroid may function specifically to enhance 
erythropoiesis and elevate oxygen transport in the blood, 
thereby helping to relieve the symptoms caused by HH 
exposure. 
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Figure 3 Regulation of WBC transcriptomic dynamics by PROG and TRIOL 

A: Principal component analysis (PCA) using genes robustly expressed in at least one of nine sample groups (i.e., no-drug treated group, PROG- 
treated group, and TRIOL-treated group in stages A, D, and F, respectively). Colors denote HH stages at which WBC samples were collected 
during HH+drug experiments, as described in Figure 1B, whereas shapes denote treatments. B: Number of differentially expressed genes (DEGs) 
between drug-treated and untreated groups in each HH stage studied. Stage A samples from both PROG and TRIOL groups were collected prior to 
drug injection; thus, DEGs of A/PROG’-A and A/TRIOL’-A reflect random variations between two groups of experimental monkeys rather than drug- 
specific effects. C: Erythrocyte-associated DEGs up-regulated by PROG. Red color and arrows denote genes up-regulated by PROG. D: 
Expression patterns and functions of TRIOL-induced DEGs in three HH-responding modules, shown as log2-transformed reads per kilobase of 
transcript per million mapped reads (RPKM) fold changes of DEGs in each sample group compared to stage A samples without TRIOL treatment. 
Light red and blue lines indicate expression patterns of individual DEGs in TRIOL-treated and untreated groups, respectively. Bold lines show mean 
expression levels of all TRIOL up- or down-regulated DEGs in each group, with error bars representing SD. UP and DOWN indicate exact number 
of up- or down-regulated DEGs by TRIOL in each module. Text below each module outlines functions associated with TRIOL-induced DEGs, with 
up and down arrows denoting gene-associated functions up- or down-regulated by TRIOL, respectively. 


Compared with PROG, TRIOL had significantly broader 
effects on the WBC transcriptome (Figure 3B), inducing 1 676 
DEGs during stages D (195 DEGs) and F (1605 DEGs). 
Interestingly, these TRIOL-induced DEGs_ overlapped 
significantly with HH-induced DEGs in all three HH-responding 
modules (Supplementary Figure S2A). Functional enrichment 
analysis revealed that TRIOL significantly up-regulated 
HIF1/hypoxia and RBC-associated inflammation, and down- 
regulated T cell receptor signaling during stage F (Figure 3D 
and Supplementary Figure S2B). Thus, in contrast to PROG, 
which had a specific function in erythropoiesis, TRIOL 
treatment altered the expression of a broader variety of genes 
in WBCs in terms of HH-related functions. 


Profound transcriptomic alterations underlie pathological 
effects of HH-induced brain injuries 
To investigate the brain gene regulatory network underlying 


the acute HH response, and the effects of PROG and TRIOL, 
we performed RNA-seq on frontal cortex samples (brains) 
from all monkeys after HH treatment (Figure 1A, B). Based on 
PCA of 14 200 robustly expressed genes, the brain samples 
could be clearly separated into the three groups: NN, HH, and 
drug treatment. The HH and NN groups could be separated at 
both the PC2 and PC3 level, as was performed for the WBC 
samples, indicating that HH was the main factor affecting the 
global brain gene expression pattern for monkeys without drug 
treatment (Figure 4A). Moreover, both agent-treated groups 
clustered tightly together, but separately from the HH group at 
the PC2 level (34.8% of total variance) and from the NN group 
at the PC3 level (23.8% of total variance) (Figure 4A), 
suggesting marked effects of both neuroprotective steroid 
agents on the brain transcriptome. 

We identified 2 992 DEGs between the NN and HH groups 
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Figure 4 Transcriptomic changes in frontal cortex in response to acute HH, PROG, and TRIOL 

A: Principal component analysis (PCA) using genes robustly expressed in at least one of four sample groups (normobaric normoxia [NN], HH, HH+ 
PROG, and HH+TRIOL). Colors denote sample group. B: Expression changes of HH-induced differentially expressed genes (DEGs) after PROG or 
TRIOL treatment, shown as log2-transformed reads per kilobase of transcript per million mapped reads (RPKM) fold changes of DEGs in each 
sample group compared to NN group. HH-induced DEGs are categorized into strongly responsive genes (SRGs), weakly responsive genes 
(WRGs), and non-responsive genes (NRGs) according to their degree of expression level recovery after drug treatment. Red and light blue lines 
denote expression patterns of individual DEGs up- and down-regulated by HH, respectively. Bold lines show mean expression levels of all HH up- 
or down-regulated DEGs in each group, with error bars representing SD. Percentages of SRGs and WRGs versus all HH-induced DEGs are 
presented above each plot. C: Expression changes of drug DEGs after PROG or TRIOL treatment, shown as log2-transformed RPKM fold changes 
of DEGs in each sample group compared to NN group. Red and light blue lines denote expression patterns of individual DEGs up- or down- 
regulated by drugs, respectively. Bold lines show mean expression levels of all DEGs up- or down-regulated by drugs in each group, with error bars 
representing SD. D: Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways enriched in HH-induced or drug-recovered DEGs. Only 
significantly enriched pathways with a false discovery rate (FDR)-adjusted P-value of <0.05 are plotted. In HH DEGs, SRGs, WRGs, and NRGs, red 
and blue dots denote pathways enriched in DEGs up- and down-regulated by HH, respectively. In drug DEGs (PROG DEGs and TRIOL DEGs), red 
and blue dots denote pathways enriched in DEGs up- and down-regulated by drugs, respectively. Color intensity of dot indicates level of 
significance, and size of dot denotes ratio of DEGs versus all expressed genes in each pathway. E: Postulated regulation of excitatory glutamate 
signaling by TRIOL. Rectangles represent proteins of TRIOL DEGs involved in glutamate signaling pathway. Blue colors denote down-regulation by 
HH. Gene symbols of DEGs associated with each protein in the pathway are listed below. F: Time-course showing intraneuronal calcium [Ca?*]; in 
primary cultured rat cortical neurons stimulated by glutamate with different doses of TRIOL. Glutamate stimulation was added at the fiftieth second 
in the experiment. 
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(see Methods), with 1446 up-regulated and 1546 down- 
regulated by HH (Supplementary Table S4), reflecting 
considerable regulatory and functional rewiring in 
brainexpressed genes in response to acute HH exposure. 
DEGs showing the largest expression fold changes (top 1% of 
all DEGs) comprised the hemoglobins (carrier of oxygen) 
HBA1 (FC=10.69, FDR=0.0000) and HBA2 (FC=4.12, 
FDR=0.000 2), and several angiogenesis-related genes, 
including adrenomedullin ADM (FC=13.58, FDR=0.000 0), 
vascular endothelial growth factor VEGFA (FC=5.48, 
FDR=0.000 0), and apelin APLN (FC=3.29, FDR=0.000 0), 
which were all up-regulated in HH brains. It is worth noting 
that many tightly associated DEGs were also significantly up- 
regulated by HH, including VEGFB (F=1.58, FDR=0.000 0), 
receptor of vascular endothelial growth factor FLT1 (FC=2.62, 
FDR=0.000 0), and receptor of adrenomedullin RAMP2 
(FC=1.45, FDR=0.006 2). This strongly suggests that an 
adaptive response to HH, by increasing local tissue 
oxygenation and increasing oxygen delivery, occurred in 
monkey brains after 52 h of HH exposure. However, VEGF 
expression can also promote BBB leakage (Schoch et al., 
2002), as observed in all monkeys exposed to HH (Figure 1F). 
Of note, CARTPT, a gene that encodes neuropeptides 
controlling feeding (Lambert et al., 1998), locomotor activity 
(Kimmel et al., 2002), and stress response (Rogge et al., 
2008), also emerged as a top DEG (FC=6.88, FDR=0.000 3) 
between NN and HH groups, and its up-regulation in HH 
brains likely played a role in inducing the symptoms of 
anorexia and motor deficits in HH monkeys. 

As with WBC transcriptome analysis, functional enrichment 
analyses revealed that HIF-1 signaling was significantly 
enriched by DEGs up-regulated in HH brains (Figure 4D), 
including EPAS1 (FC=1.37, FDR=0.0285) and EGLN1 
(FC=1.48, FDR=0.004 7), two genes strongly associated with 
high-altitude adaptation in highland populations (Simonson et 
al., 2010; Yi et al., 2010), and H/F3A (FC=4.36, FDR=0.000 0), 
an important regulator of HIF1A and EPAS7 involved in 
adaptive responses to hypoxia (Ravenna et al., 2015). 
Moreover, given that the top DEGs (e.g., ADM, VEGF, APLN, 
HBA1, and HBA2) are regulated by HIF-1 signaling (Chen et 
al., 2012; Eyries et al., 2008; Ramakrishnan et al., 2014; Saha 
et al., 2014), this suggests that HIF-1 signaling also plays a 
main regulatory role in adaptive responses to acute HH in 
brains. Additionally, genes involved in the complement and 
coagulation cascades, and in leukocyte transendothelial 
migration, were up-regulated in HH-exposed brains (Figure 4D 
and Supplementary Figure S3A) and WBCs (Figure 2D), 
implying a link between peripheral inflammation and 
neuroinflammation in the CNS. 

As a consequence of the decrease in oxygen supply under 
HH conditions, several pathways involved in aerobic 
respiration were found to be down-regulated in HH brains, 
including oxidative phosphorylation, respiratory electron 
transport chain, tricarboxylic acid (TCA) cycle, and pyruvate 
metabolism. On the other hand, pathways involved in 
anaerobic respiration were up-regulated, including glycolysis 


and mTOR signaling (Figure 4D and Supplementary Figure 
S3A). These findings suggest that, regardless of the adaptive 
response by stimulating angiogenesis and increasing oxygen 
delivery, the monkey brains suffered an energy deficiency, 
which may have induced autophagy and neuronal injuries 
through mTOR signaling (Dunlop & Tee, 2014). 

Protein misfolding and endoplasmic reticulum (ER) stress 
are involved in many types of brain injury and 
neurodegenerative diseases (Giffard et al., 2004; Soto & 
Estrada, 2008). Indeed, we found that many key genes 
involved in the protein folding process were down-regulated by 
HH (Supplementary Figure S3A), including those encoding 
members of the HSP40 family, prefoldins, and FK506 binding 
proteins (Supplementary Table S4). Members of the HSP40 
family are crucial partners for HSP70 chaperones (Qiu et al., 
2006), and are important for hypoxic tolerance (Jain et al., 
2014) and injury resistance (Hasegawa et al., 2018) in the 
CNS. Thus, their down-regulation may exacerbate HH- 
induced brain injury. 

Taken together, our results revealed profound neural gene 
expression alterations after exposure to acute HH, and 
rewiring of genetic pathways, which may contribute to an 
adaptive response to acute HH and underlie the pathological 
effects of the associated brain injuries. 


Brain transcriptomic changes induced by HH are reversed 
by treatment with PROG and TRIOL 

We then assessed how PROG and TRIOL affect the 
expression of HH-altered genes (see Methods). Of all HH- 
altered genes, the expression levels were restored to normal 
(defined by reference to the NN group) in 50.67% and 48.09% 
of cases by PROG and TRIOL treatment, respectively, i. e., 
‘strongly responsive genes' (SRGs) (Figure 4B). Furthermore, 
31.52% and 29.45% of HH-altered genes were partially 
restored by PROG and TRIOL treatment, respectively, i. e., 
‘weakly responsive genes' (WRGs) (Figure 4B). In addition, 
the expression levels of 17.15% and 20.82% of HH-altered 
genes were unaffected by PROG and TRIOL treatment, 
respectively, i.e., 'non-responsive genes' (NRGs) (Figure 4B). 
These results indicate that the abnormal expression levels of 
most HH-induced DEGs were completely or partially restored 
by PROG (82.19%) or TRIOL (77.54%) treatment. It should be 
noted that the SRGs and WRGs associated with the two drug 
treatments largely overlapped with each other (Supplementary 
Figure S3B), suggesting that PROG and TRIOL have broadly 
similar effects on reversing WHH-induced transcriptomic 
changes in the brain. 

PROG and TRIOL SRGs were particularly enriched in 
components of the key pathways of aerobic respiration (Figure 
4D and Supplementary Figure S3A), suggesting an 
improvement in aerobic energy supply in response to PROG 
and TRIOL treatment. Specifically, of the 44 genes involved in 
oxidative phosphorylation and down-regulated by HH, 25 
(56.8%) were PROG SRGs and 39 (88.6%) were TRIOL 
SRGs (Fisher's exact test, P=0.002). This implies that TRIOL 
may be more effective at improving the aerobic energy supply 
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than PROG. The abnormal HH-induced expression of genes 
with roles in protein folding was also restored by both PROG 
and TRIOL treatment (Figure 4D and Supplementary Figure 
S3A), which may relieve neuronal injuries induced by protein 
folding and ER stress. In addition, the up-regulation of the 
behavior-related gene CARTPT by HH was also attenuated by 
both drugs to a certain degree (Supplementary Table S4), 
consistent with the alleviation of anorexia and motor deficits in 
drug-treated monkeys. 

Of note, the up-regulations of several top DEGs involved in 
angiogenesis by HH were not recovered by PROG or TRIOL, 
including ADM, VWF, and VEGFA, suggesting that 
angiogenesis was also active in drug-treated groups 
(Supplementary Table S4). Intriguingly, PROG but not TRIOL 
further up-regulated several key genes involved in 
angiogenesis and oxygen delivery, including VEGFA (HH-vs- 
PROG: FC=1.23, FDR=0.032 3), its receptor FLT1 (HH-vs- 
PROG: FC=1.23, FDR=0.032 3), and hemoglobins HBB (HH- 
vs-PROG: FC=2.33, FDR=0.002 4) and HBA2 (HH-vs-PROG: 
FC=1.97, FDR=0.000 0) (Supplementary Table S4). These 
results further emphasize the potential role of PROG in 
relieving HH-induced symptoms by elevating oxygen transport 
in both WBCs and brains. 


Gene regulatory mechanisms by which TRIOL suppresses 
glutamate-mediated excitotoxicity 
In the brain transcriptome, we detected many genes that were 
specifically rewired by either PROG or TRIOL, even though 
their expression levels were not altered by HH treatment (see 
Methods, Figure 4C). These included 380 and 743 genes 
whose expression levels were induced by PROG and TRIOL, 
respectively. Intriguingly, among them, we found that TRIOL, 
but not PROG, significantly induced the enrichment of DEGs 
in many key excitatory neuronal signaling pathways (Figure 
4D and Supplementary Figure S3A), including the 
glutamatergic synapse, L-glutamate import, glutamate 
receptor, and calcium signaling pathways. This suggests that 
TRIOL may have a specific effect on the regulation of 
glutamate-induced excitotoxicity. Glutamate induces 
excitotoxicity by promoting the excessive influx of intracellular 
calcium [Ca?*] through ionotropic glutamate receptors, 
particularly through the N-methyl-d-aspartate (NMDA) 
receptor. This is the main cause of hypoxia-induced neuronal 
death in many hypoxia-induced brain injuries (Wroge et al., 
2012). Many key genes involved in this process were down- 
regulated by TRIOL, including NMDA receptors GRIN2A 
(FC=1.79, FDR=0.000 0)and GRIN2B(FC=2.23, FDR=0.000 3), 
vesicular glutamate transporter (VGLUT) SLC17A7 (FC=1.63, 
FDR=0.000 0), excitatory amino acid transporter (EAAT) 
SLC1A2 (FC=1.60, FDR=0.018 4), group | metabotropic 
glutamate receptor GRM5 (FC=1.48, FDR=0.030 4), and its 
downstream adenylate cyclase family members ADCY1 
(FC=1.52, FDR=0.029 1), ADCY3 (FC=1.39, FDR 0.009 8), 
ADCY5 (FC=1.55, FDR 0.003 8), and ADCY9 (FC=1.59, 
FDR=0.000 2) (Figure 4E and Supplementary Table S4). 
Glutamate-induced intracellular calcium ([Ca?*],) overload is 
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part of the signature indicating excitotoxicity in neurons 
(Scannevin & Huganir, 2000). To further examine whether 
TRIOL can indeed reduce neuronal excitotoxicity by directly 
suppressing the intracellular accumulation of calcium, we 
measured glutamate-induced [Ca?*]i in primary rat cortical 
neurons with and without TRIOL treatment. We observed that 
stimulation with 100 pmol/L glutamate alone dramatically 
increased [Ca] in neurons within 50s, whereas 
supplementation with TRIOL decreased these glutamate- 
induced effects in a dose-dependent fashion (Figure 4F). 
Taken together, our results suggest that TRIOL may suppress 
excitotoxicity via the indirect transcriptional regulation of key 
glutamate signaling-related genes through direct modulation of 
membrane permeability for calcium. 


DISCUSSION 


In this study, we demonstrated that cynomolgus monkeys 
exposed to acute HH exhibited various symptoms, including 
anorexia, vomiting, motor deficits, and ataxia, as observed in 
human patients suffering acute high-altitude diseases and 
which are mostly attributed to pathological changes in the 
brain (Wilson et al., 2009). The pathophysiological diagnoses 
in the brain suggest that cynomolgus monkeys developed 
severe cerebral edema, similar to that found in human HACE 
patients (Hackett & Roach, 2004). Hence, these data suggest 
cynomolgus monkeys can be used as non-human primate 
models for studying the disorders of oxygen homeostasis. 
Moreover, the genetic regulatory networks shown in this study 
to be involved in the development of hypoxia-induced brain 
malfunctions provide valuable insights into the relationship 
between the mechanism of oxygen homeostasis regulation 
and related diseases. Our study confirmed the fundamental 
role of several classical pathways, such as HIF-1 signaling, in 
the organism level response to hypoxia. Moreover, we also 
revealed that several molecules, such as VDR, serve as 
central connectors of different acute HH response pathways; 
these have been largely ignored in previous studies using 
human cell lines or mice as experimental models. 

The dynamic transcriptome of WBCs from non-human 
primates may help to improve our understanding of the 
progression of HACE. First of all, most genes were 
categorized as belonging to the impulse module, including 
those involved in HIF-1 signaling; the expression of these 
genes was stimulated during early HH stages but recovered in 
late stages, suggesting that hypoxia is indeed the primary 
cause of injury under acute HH conditions. The hub genes in 
the impulse module may account for the subsequent 
expression changes in genes in the late stages of HH and 
may be responsible for the fundamental pathophysiological 
changes involved in acclimatization or injury. For example, 
RBC-associated inflammation is a novel pathway in late 
stages and may interact with the progress of HACE. Our 
transcriptomic analyses combining WBCs and brains also 
highlighted a complex organismal response to acute HH, with 
active crosstalk between the acclimatization that maintains 


oxygen homeostasis and adverse pathological feedback 
caused by inflammation. Specifically, in addition to the 
activation of the HIF-1 signaling pathway, consistent with 
previous reports (Semenza, 2009; Wilson et al., 2009), we 
also found that HH may induce a hierarchy of transcriptional 
programs of both innate and adaptive immunity, while 
progressively promoting inflammation in both the blood and 
brain. 

The development of treatments for brain injury associated 
with HH, such as HACE, remains a formidable challenge. Our 
pharmacogenomic studies reported here provide preliminary 
but nevertheless valuable insights into the feasibility of a 
generalized neuroprotection strategy. Further, they present 
possible drug candidates, PROG and TRIOL, which both 
alleviated acute HH-induced brain injuries in our primate 
model. PROG is reported to be neuroprotective by stimulating 
breathing, increasing oxygen-carrying capacity, and 
modulating the GABAa receptor to suppress excitotoxicity 
(Singh & Su, 2013), whereas our results suggested a role in 
stimulating erythropoiesis and increasing oxygen delivery. 
Medroxyprogesterone (MPA), a synthetic analogue of PROG 
without neuroprotectant properties (Singh & Su, 2013), is used 
to treat chronic mountain sickness by stimulating breathing 
and increasing blood oxygen levels (Wright et al., 2004). 
However, MPA has no significant prophylactic effect for acute 
mountain sickness (AMS) (Wright et al., 2004), which further 
emphasizes the importance of neuroprotection in treating 
HACE. It is intriguing that the two neuroprotectants, PROG 
and TRIOL, appear to have different modes of action and act 
through distinct functional mechanisms to rewire the 
expression of genes in the blood and brain. Our previous 
study reported that TRIOL can protect primary cortical 
neurons from hypoxia/reoxygenation injury by maintaining 
mitochondrial functions (Chen et al., 2013a). Here, we further 
revealed that the neuroprotective effect of TRIOL may be 
mediated by attenuation of glutamate-induced excitotoxicity. 

In conclusion, our research extends the transcriptomic 
analysis of acute HH responses at an integrated spatial and 
temporal level in non-human primates. We presented a kinetic 
transcriptomic analysis of WBCs and the transcriptomic profile 
of brains after exposure to acute HH. Our data revealed 
thousands of genes regulated by acute HH and categorized to 
different co-expression modules depending on their spatial 
and temporal dimensions. Our study also revealed that 
neuroprotective steroids such as PROG and TRIOL have 
significant recovery effects on acute HH-induced behavioral 
deficits and brain damage, suggesting the potential of using 
novel neuroprotectants as therapeutic drugs for acute high- 
altitude diseases. However, we also acknowledge that there 
are ~25 million years of evolutionary divergence between 
cynomolgus monkeys and humans (Yan et al., 2011), and 
some routine physiological data such as arterial blood 
pressure, arterial oxygenation, and ventilation were not 
collected in this study due to the limitation of experimental 
conditions. Thus, we cannot rule out that the pathological and 
gene regulatory changes observed in the cynomolgus 


monkeys may differ from that in humans suffering HH-induced 
disorders. Future studies with relevant physiological 
measurements could help better understand how animals 
respond to severe hypoxia, and how neuroprotective steroids 
function to protect HH-induced brain damage. 
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